feat: add vaf calculation for strelka results #109
Conversation
WalkthroughThis pull request introduces a new Conda environment configuration by adding a YAML file that specifies the channels and a dependency on Changes
Sequence Diagram(s)sequenceDiagram
participant W as Workflow Engine
participant R as calculate_vaf Rule
participant E as Conda Environment (cyvcf.yaml)
participant S as calc-vaf.py Script
W->>R: Trigger calculate_vaf rule with input BCF file
R->>E: Activate Conda environment
R->>S: Execute calc-vaf.py script
S->>S: Call get_snv_allele_freq / get_indel_allele_freq
S-->>R: Return calculated VAF and output BCF file
R-->>W: Provide results and logs
sequenceDiagram
participant C as Caller
participant F as get_vaf_status Function
C->>F: Invoke get_vaf_status with wildcards (including vaf_benchmark)
alt vaf_benchmark == "tbc"
F-->>C: Return True immediately
else vaf_benchmark != "tbc"
F-->>C: Evaluate existing VAF checks
end
Assessment against linked issues
Suggested reviewers
Poem
Note ⚡️ AI Code Reviews for VS Code, Cursor, WindsurfCodeRabbit now has a plugin for VS Code, Cursor and Windsurf. This brings AI code reviews directly in the code editor. Each commit is reviewed immediately, finding bugs before the PR is raised. Seamless context handoff to your AI code agent ensures that you can easily incorporate review feedback. Note ⚡️ Faster reviews with cachingCodeRabbit now supports caching for code and dependencies, helping speed up reviews. This means quicker feedback, reduced wait times, and a smoother review experience overall. Cached data is encrypted and stored securely. This feature will be automatically enabled for all accounts on May 16th. To opt out, configure ✨ Finishing Touches
Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. 🪧 TipsChatThere are 3 ways to chat with CodeRabbit:
SupportNeed help? Create a ticket on our support page for assistance with any issues or questions. Note: Be mindful of the bot's finite context window. It's strongly recommended to break down tasks such as reading entire modules into smaller chunks. For a focused discussion, use review comments to chat about specific files and their changes, instead of using the PR comments. CodeRabbit Commands (Invoked using PR comments)
Other keywords and placeholders
CodeRabbit Configuration File (
|
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 2
🧹 Nitpick comments (4)
workflow/scripts/calc-vaf.py (1)
3-4: Remove hardcoded file paths from comments.The commented file paths appear to be local development paths. Remove them to avoid confusion.
-#vcf = "/Users/famke/01-pm4onco/osf-download/pipeline-results-of-imgag-data/qbic/strelka/tumor_5perc_vs_normal_5perc.strelka.somatic_snvs_VEP.ann.vcf" #snakemake.input -#indel = "/Users/famke/01-pm4onco/osf-download/pipeline-results-of-imgag-data/qbic/strelka/tumor_5perc_vs_normal_5perc.strelka.somatic_indels_VEP.ann.vcf.gz"workflow/rules/eval.smk (1)
94-105: Enhance the calculate_vaf rule configuration.The rule needs some improvements:
- Add resources section for memory management
- Fix the log path to include {wildcards.callset}
- Add benchmark section for performance tracking
Here's a suggested implementation:
rule calculate_vaf: input: "results/filtered-variants/{callset}.bcf" output: "results/calculate-vaf/{callset}.added-vaf.bcf" log: - "logs/calculate-vaf/" + "logs/calculate-vaf/{callset}.log" + benchmark: + "benchmarks/calculate-vaf/{callset}.tsv" + resources: + mem_mb=4000 conda: "../envs/cyvcf.yaml" script: "../scripts/calc-vaf.py"workflow/rules/common.smk (1)
458-459: Document the "tbc" condition in get_vaf_status.Add a comment explaining what "tbc" means and why it triggers a True return value.
if vaf_benchmark is None: return False + # Return True for "tbc" (to be calculated) to enable VAF calculation + # for variants that don't have pre-calculated VAF values if vaf_benchmark == "tbc": return Trueworkflow/envs/cyvcf.yaml (1)
1-5: Enhance the conda environment configuration.The configuration needs improvements:
- Pin the cyvcf2 version for reproducibility
- Add a newline at the end of file
Here's a suggested implementation:
channels: - conda-forge - bioconda dependencies: - - cyvcf2 + - cyvcf2=0.30.22 +🧰 Tools
🪛 yamllint (1.35.1)
[error] 5-5: no new line character at the end of file
(new-line-at-end-of-file)
📜 Review details
Configuration used: CodeRabbit UI
Review profile: CHILL
Plan: Pro
📒 Files selected for processing (5)
.gitignore(1 hunks)workflow/envs/cyvcf.yaml(1 hunks)workflow/rules/common.smk(1 hunks)workflow/rules/eval.smk(1 hunks)workflow/scripts/calc-vaf.py(1 hunks)
✅ Files skipped from review due to trivial changes (1)
- .gitignore
🧰 Additional context used
🪛 yamllint (1.35.1)
workflow/envs/cyvcf.yaml
[error] 5-5: no new line character at the end of file
(new-line-at-end-of-file)
workflow/scripts/calc-vaf.py
Outdated
| def get_indel_allele_freq(vcf): | ||
| for variant in VCF(vcf): | ||
| tier1RefCounts = variant.format("TAR")[0,0] | ||
| tier1AltCounts = variant.format("TIR")[0,0] | ||
|
|
||
| vaf = tier1AltCounts / (tier1AltCounts + tier1RefCounts) | ||
|
|
||
| print(vaf) |
There was a problem hiding this comment.
🛠️ Refactor suggestion
Apply similar improvements to INDEL VAF calculation.
The INDEL function needs the same enhancements as the SNV function.
Here's a suggested implementation:
def get_indel_allele_freq(vcf):
+ """Calculate Variant Allele Frequency (VAF) for INDELs.
+
+ Args:
+ vcf: Path to the VCF file containing INDEL variants
+
+ Returns:
+ List of tuples containing (variant_id, vaf)
+ """
+ results = []
for variant in VCF(vcf):
tier1RefCounts = variant.format("TAR")[0,0]
tier1AltCounts = variant.format("TIR")[0,0]
- vaf = tier1AltCounts / (tier1AltCounts + tier1RefCounts)
+ total_counts = tier1AltCounts + tier1RefCounts
+ if total_counts == 0:
+ vaf = 0.0
+ else:
+ vaf = tier1AltCounts / total_counts
- print(vaf)
+ results.append((variant.ID or f"{variant.CHROM}:{variant.POS}", vaf))
+ return results📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| def get_indel_allele_freq(vcf): | |
| for variant in VCF(vcf): | |
| tier1RefCounts = variant.format("TAR")[0,0] | |
| tier1AltCounts = variant.format("TIR")[0,0] | |
| vaf = tier1AltCounts / (tier1AltCounts + tier1RefCounts) | |
| print(vaf) | |
| def get_indel_allele_freq(vcf): | |
| """Calculate Variant Allele Frequency (VAF) for INDELs. | |
| Args: | |
| vcf: Path to the VCF file containing INDEL variants | |
| Returns: | |
| List of tuples containing (variant_id, vaf) | |
| """ | |
| results = [] | |
| for variant in VCF(vcf): | |
| tier1RefCounts = variant.format("TAR")[0,0] | |
| tier1AltCounts = variant.format("TIR")[0,0] | |
| total_counts = tier1AltCounts + tier1RefCounts | |
| if total_counts == 0: | |
| vaf = 0.0 | |
| else: | |
| vaf = tier1AltCounts / total_counts | |
| results.append((variant.ID or f"{variant.CHROM}:{variant.POS}", vaf)) | |
| return results |
workflow/scripts/calc-vaf.py
Outdated
| def get_snv_allele_freq(vcf): | ||
| for variant in VCF(vcf): | ||
| refCounts = variant.format(variant.REF + "U") | ||
| altCounts = variant.format(variant.ALT[0] + "U") | ||
|
|
||
| # TODO: check which value is the correct one from the matrix (this leads to many zero VAF) | ||
| tier1RefCounts = refCounts[0, 0] | ||
| tier1AltCounts = altCounts[0, 0] | ||
|
|
||
| vaf = tier1AltCounts / (tier1AltCounts + tier1RefCounts) | ||
|
|
||
| print(vaf) |
There was a problem hiding this comment.
🛠️ Refactor suggestion
Enhance the SNV VAF calculation implementation.
The current implementation needs improvements:
- Add error handling for division by zero
- Add docstring explaining the function
- Store results instead of just printing
- Validate the matrix values as noted in TODO
Here's a suggested implementation:
def get_snv_allele_freq(vcf):
+ """Calculate Variant Allele Frequency (VAF) for SNVs.
+
+ Args:
+ vcf: Path to the VCF file containing SNV variants
+
+ Returns:
+ List of tuples containing (variant_id, vaf)
+ """
+ results = []
for variant in VCF(vcf):
refCounts = variant.format(variant.REF + "U")
altCounts = variant.format(variant.ALT[0] + "U")
- # TODO: check which value is the correct one from the matrix (this leads to many zero VAF)
tier1RefCounts = refCounts[0, 0]
tier1AltCounts = altCounts[0, 0]
- vaf = tier1AltCounts / (tier1AltCounts + tier1RefCounts)
+ total_counts = tier1AltCounts + tier1RefCounts
+ if total_counts == 0:
+ vaf = 0.0
+ else:
+ vaf = tier1AltCounts / total_counts
- print(vaf)
+ results.append((variant.ID or f"{variant.CHROM}:{variant.POS}", vaf))
+ return results📝 Committable suggestion
‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.
| def get_snv_allele_freq(vcf): | |
| for variant in VCF(vcf): | |
| refCounts = variant.format(variant.REF + "U") | |
| altCounts = variant.format(variant.ALT[0] + "U") | |
| # TODO: check which value is the correct one from the matrix (this leads to many zero VAF) | |
| tier1RefCounts = refCounts[0, 0] | |
| tier1AltCounts = altCounts[0, 0] | |
| vaf = tier1AltCounts / (tier1AltCounts + tier1RefCounts) | |
| print(vaf) | |
| def get_snv_allele_freq(vcf): | |
| """Calculate Variant Allele Frequency (VAF) for SNVs. | |
| Args: | |
| vcf: Path to the VCF file containing SNV variants | |
| Returns: | |
| List of tuples containing (variant_id, vaf) | |
| """ | |
| results = [] | |
| for variant in VCF(vcf): | |
| refCounts = variant.format(variant.REF + "U") | |
| altCounts = variant.format(variant.ALT[0] + "U") | |
| tier1RefCounts = refCounts[0, 0] | |
| tier1AltCounts = altCounts[0, 0] | |
| total_counts = tier1AltCounts + tier1RefCounts | |
| if total_counts == 0: | |
| vaf = 0.0 | |
| else: | |
| vaf = tier1AltCounts / total_counts | |
| results.append((variant.ID or f"{variant.CHROM}:{variant.POS}", vaf)) | |
| return results |
There was a problem hiding this comment.
Actionable comments posted: 1
♻️ Duplicate comments (2)
workflow/scripts/calc-vaf.py (2)
6-16: Add error handling for division by zero in SNV VAF calculation.The function doesn't handle the case where both reference and alternate counts are zero, which would cause a division by zero error. Additionally, a docstring would improve code readability.
def get_snv_allele_freq(variant): + """Calculate Variant Allele Frequency (VAF) for SNVs. + + Args: + variant: A cyvcf2 variant object containing SNV data + + Returns: + float: The calculated VAF + """ refCounts = variant.format(variant.REF + "U") altCounts = variant.format(variant.ALT[0] + "U") # TODO: check which value is the correct one from the matrix (this leads to many zero VAF) tier1RefCounts = refCounts[0, 0] tier1AltCounts = altCounts[0, 0] - vaf = tier1AltCounts / (tier1AltCounts + tier1RefCounts) + total_counts = tier1AltCounts + tier1RefCounts + if total_counts == 0: + vaf = 0.0 + else: + vaf = tier1AltCounts / total_counts return vaf
18-23: Add error handling for division by zero in INDEL VAF calculation.Similar to the SNV function, this function lacks error handling for division by zero and would benefit from a proper docstring.
def get_indel_allele_freq(variant): + """Calculate Variant Allele Frequency (VAF) for INDELs. + + Args: + variant: A cyvcf2 variant object containing INDEL data + + Returns: + float: The calculated VAF + """ tier1RefCounts = variant.format("TAR")[0,0] tier1AltCounts = variant.format("TIR")[0,0] - vaf = tier1AltCounts / (tier1AltCounts + tier1RefCounts) + total_counts = tier1AltCounts + tier1RefCounts + if total_counts == 0: + vaf = 0.0 + else: + vaf = tier1AltCounts / total_counts return vaf
🧹 Nitpick comments (1)
workflow/scripts/calc-vaf.py (1)
1-4: Remove unused importargparse.The
argparsemodule is imported but not used anywhere in the script. Since there's no command line argument parsing functionality implemented, this import should be removed.from cyvcf2 import VCF, Writer import sys -import argparse import numpy as np🧰 Tools
🪛 Ruff (0.8.2)
3-3:
argparseimported but unusedRemove unused import:
argparse(F401)
📜 Review details
Configuration used: CodeRabbit UI
Review profile: CHILL
Plan: Pro
Cache: Disabled due to data retention organization setting
Knowledge Base: Disabled due to data retention organization setting
📒 Files selected for processing (1)
workflow/scripts/calc-vaf.py(1 hunks)
🧰 Additional context used
🪛 Ruff (0.8.2)
workflow/scripts/calc-vaf.py
3-3: argparse imported but unused
Remove unused import: argparse
(F401)
⏰ Context from checks skipped due to timeout of 90000ms (1)
- GitHub Check: Testing
🔇 Additional comments (3)
workflow/scripts/calc-vaf.py (3)
25-78: LGTM! Robust VAF calculation function with proper error handling.The
calculate_vaffunction is well implemented with:
- Comprehensive docstring explaining inputs and outputs
- Proper error handling for missing AD field
- NaN values for missing or zero-divisor cases
- Efficient numpy array operations
- Clear comments explaining the code logic
80-144: LGTM! Well-structured main function with thorough error handling.The
add_vaf_to_vcffunction is well-designed with:
- Good error handling for file operations and format field setting
- Robust processing that continues even if some variants have issues
- Progress logging for large files
- Clear structure and documentation
6-23: Clarify usage of unused variant-specific VAF functions.The
get_snv_allele_freqandget_indel_allele_freqfunctions are defined but not used in the script. The more genericcalculate_vafis used instead. Consider either:
- Removing these unused functions if they're not needed
- Adding a comment explaining their purpose if they're kept for future use or reference
- Integrating them into the main workflow if they're meant to be used for Strelka-specific calculations
Could you clarify the intended purpose of these two functions in relation to the overall VAF calculation workflow?
| print(f"Warning: AD field missing for variant at {variant.CHROM}:{variant.POS}. Skipping VAF calculation.", file=sys.stderr) | ||
| return None | ||
|
|
||
| n_samples = len(samples) | ||
| n_alleles = len(variant.alleles) # Includes reference | ||
| n_alt_alleles = n_alleles - 1 | ||
|
|
||
| # Initialize VAF array with NaNs | ||
| # Shape: (n_samples, n_alt_alleles) | ||
| vaf_values = np.full((n_samples, n_alt_alleles), np.nan, dtype=np.float32) | ||
|
|
||
| for i in range(n_samples): | ||
| sample_ad = ad[i] | ||
|
|
||
| # Check for missing AD data for the sample (represented by negative values or could be None depending on VCF) | ||
| # cyvcf2 often uses negative numbers for missing integers in FORMAT fields like AD | ||
| if np.any(sample_ad < 0): | ||
| # Keep VAF as NaN if AD is missing for this sample | ||
| continue | ||
|
|
||
| # Calculate total depth for this sample | ||
| total_depth = np.sum(sample_ad) | ||
|
|
||
| if total_depth == 0: | ||
| # Avoid division by zero, keep VAFs as NaN | ||
| continue | ||
|
|
||
| # Calculate VAF for each alternate allele | ||
| for j in range(n_alt_alleles): | ||
| alt_depth = sample_ad[j + 1] # AD format is [ref_depth, alt1_depth, alt2_depth, ...] | ||
| vaf = alt_depth / total_depth | ||
| vaf_values[i, j] = vaf | ||
|
|
||
| return vaf_values | ||
|
|
||
| def add_vaf_to_vcf(input_vcf_path, output_vcf_path): | ||
| """ | ||
| Reads an input VCF, calculates VAF for each variant/sample, adds it | ||
| as a new FORMAT field 'VAF', and writes to an output VCF file. | ||
| """ | ||
| # Open the input VCF file | ||
| try: | ||
| vcf_reader = VCF(input_vcf_path) | ||
| except Exception as e: | ||
| print(f"Error opening input VCF file '{input_vcf_path}': {e}", file=sys.stderr) | ||
| sys.exit(1) | ||
|
|
||
|
|
||
| # Add the new VAF FORMAT field definition to the header | ||
| # Number='A' means one value per alternate allele | ||
| # Type='Float' for the frequency value | ||
| try: | ||
| vcf_reader.add_format_to_header({ | ||
| 'ID': 'VAF', | ||
| 'Description': 'Variant Allele Frequency calculated from AD field (Alt Depth / Total Depth)', | ||
| 'Type': 'Float', | ||
| 'Number': 'A' # One value per alternate allele | ||
| }) | ||
| except ValueError as e: | ||
| # Catch error if the field already exists | ||
| print(f"Warning: FORMAT field 'VAF' might already exist in header: {e}", file=sys.stderr) | ||
|
|
||
|
|
||
| # Create a VCF writer object using the modified header | ||
| try: | ||
| vcf_writer = Writer(output_vcf_path, vcf_reader) | ||
| except Exception as e: | ||
| print(f"Error creating output VCF file '{output_vcf_path}': {e}", file=sys.stderr) | ||
| vcf_reader.close() | ||
| sys.exit(1) | ||
|
|
||
| print(f"Processing VCF: {input_vcf_path}") | ||
| print(f"Writing output to: {output_vcf_path}") | ||
|
|
||
| processed_count = 0 | ||
| # Iterate through each variant in the VCF | ||
| for variant in vcf_reader: | ||
| # Calculate VAFs for all samples for the current variant | ||
| vaf_array = calculate_vaf(variant, vcf_reader.samples) | ||
|
|
||
| # Add the calculated VAFs to the variant's FORMAT fields | ||
| # The vaf_array must have shape (n_samples, n_alt_alleles) | ||
| if vaf_array is not None: | ||
| try: | ||
| # Use set_format to add/update the VAF field for all samples | ||
| variant.set_format('VAF', vaf_array) | ||
| except Exception as e: | ||
| print(f"Error setting VAF format for variant at {variant.CHROM}:{variant.POS}: {e}", file=sys.stderr) | ||
| # Decide if you want to skip writing this variant or write without VAF | ||
| # Here, we'll still write the variant but VAF might be missing/incorrect | ||
|
|
||
| # Write the (potentially modified) variant record to the output file | ||
| vcf_writer.write_record(variant) | ||
| processed_count += 1 | ||
| if processed_count % 1000 == 0: | ||
| print(f"Processed {processed_count} variants...", file=sys.stderr) | ||
|
|
||
| # Close the VCF reader and writer | ||
| vcf_reader.close() | ||
| vcf_writer.close() |
There was a problem hiding this comment.
🛠️ Refactor suggestion
Add a main block to make the script executable.
The script doesn't have a way to be executed from the command line. Since the file will likely be used as a command-line tool in a Snakemake workflow, you should add a main block to parse arguments and execute the add_vaf_to_vcf function.
import sys
import numpy as np
+import argparse
...existing code...
vcf_reader.close()
vcf_writer.close()
print(f"Finished processing. Total variants processed: {processed_count}")
+if __name__ == "__main__":
+ parser = argparse.ArgumentParser(description="Add VAF information to VCF files")
+ parser.add_argument("input_vcf", help="Input VCF file")
+ parser.add_argument("output_vcf", help="Output VCF file with VAF added")
+ args = parser.parse_args()
+
+ add_vaf_to_vcf(args.input_vcf, args.output_vcf)Committable suggestion skipped: line range outside the PR's diff.
🧰 Tools
🪛 Ruff (0.8.2)
3-3: argparse imported but unused
Remove unused import: argparse
(F401)
🤖 Prompt for AI Agents
In workflow/scripts/calc-vaf.py around lines 1 to 144, the script lacks a main
block to allow command-line execution. Add a main block at the end of the file
that uses argparse to parse input and output VCF file paths from command-line
arguments, then calls the add_vaf_to_vcf function with those arguments. This
will enable the script to be run directly and integrated into workflows like
Snakemake.
Closes #100
This is a work in progress. The following steps still have to be solved:
Summary by CodeRabbit